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Abstract. To understand the nonlinear dynamics of the Parker scenario for coronal heating, long- 
time high-resolution simulations of the dynamics of a coronal loop in cartesian geometry are carried 
out. A loop is modeled as a box extended along the direction of the strong magnetic field Bq in which 
the system is embedded. At the top and bottom plates, which represent the photosphere, velocity 
fields mimicking photospheric motions are imposed. 

We show that the nonlinear dynamics is described by different regimes of MHD anisotropic 
turbulence, with spectra characterized by intertial range power laws whose indexes range from 
Kolmogorov-like values (~ 5/3) up to ~ 3. We briefly describe the bearing for coronal heating 
rates. 
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INTRODUCTION 

Coronal heating is one of the outstanding problems in solar physics. Although the 
correlation of coronal activity with the intensity of photospheric magnetic fields seems 
beyond doubt, and there is large agreement that photospheric motions are the source of 
the energy flux that sustains an active region (~ 10 7 ergcm~ 2 s~ l ), the debate currently 
focuses on the physical mechanisms responsible for the transport, storage and dissipation 
(i.e. conversion to heat and/or particle acceleration) of this energy from the photosphere 
to the corona. 

A promising model is that proposed by Parker [1,2], who was the first to suggest that 
coronal heating could be the necessary outcome of an energy flux associated with the 
tangling of coronal field lines by photospheric motions. 

Over the years a number of numerical experiments have been carried out to investigate 
the Parker problem, with particular emphasis on exploring how the shuffling of magnetic 
field line footpoints leads to current sheet formation, and to estimate the heating rate. 

3D cartesian simulations have been performed by Mikic et al. [3], Longcope and 
Sudan [4], Hendrix and Van Hoven [5], and Dmitruk and Gomez [6]. A complex coronal 
magnetic field results from the photospheric field line random walk, and though the field 
does not, strictly speaking, evolve through a sequence of static force-free equilibrium 
states (the original Parker hypothesis), magnetic energy nonetheless tends to dominate 
kinetic energy in the system. In this limit the field is structured by current sheets 
elongated along the axial direction. The results from these studies agreed qualitatively 
among themselves, in that all simulations display the development of field aligned 



current sheets. However, estimates of the dissipated power and its scaling characteristics 
differed largely, depending on the way in which extrapolations from low to large values 
of the plasma conductivity of the properties such as inertial range power law indices 
were carried out. 

The low resolution of the previous 3D studies has been partially overcome by 2D 
numerical simulations of incompressible MHD with magnetic forcing (Einaudi et al. 
[7], Georgoulis et al. [8], Dmitruk et al. [9], Einaudi and Velli [10]), which showed that 
turbulent current sheets dissipation is distributed intermittently, and that the statistics 
of dissipation events, in terms of total energy, peak energy and event duration displays 
power laws not unlike the distribution of observed emission events in optical, ultraviolet 
and x-ray wavelengths of the quiet solar corona. 

More recently a first attempt to simulate full 3D sections of the solar corona with a 
realistic geometry has been performed by Gudiksen and Nordlund [11]. At the moment 
the very low resolution attainable with this kind of simulations does not allow the 
development of turbulence. The transfer of energy from the scale of convection cells 
~ 1000 km toward smaller scales is in fact inhibited, because the smaller scales are not 
resolved (their linear resolution is in fact ~ 500 km). 

While in the future these global simulations will be able to reach the necessary high 
resolutions, to investigate the nonlinear dynamics of the Parker scenario at relatively 
high Reynolds numbers, we have recently performed high-resolution long-time simula- 
tion of the aforementioned cartesian model (Rappazzo et al. [12]). 

In the next sections we describe the coronal loop model, the simulations we have 
carried out, and give simple scaling arguments to understand the energy spectral slopes. 

PHYSICAL MODEL 

A coronal loop is a closed magnetic structure threaded by a strong axial field, with the 
footpoints rooted in the photosphere. This makes it a strongly anisotropic system, as 
measured by the relative magnitude of the Alfven velocity associated with the axial 
magnetic field va ~ 2000 kms -1 compared to the typical photospheric velocity u p h ~ 
1 kms -1 . This means that the relative amplitude of the Alfven waves that are launched 
into the corona is very small and, as an efficient energy cascade takes place [12], the 
relative amplitude of the fields which develop in the orthogonal planes remains small 
compared to the dominant axial magnetic field. 

We study the loop dynamics in a simplified cartesian geometry, neglecting any cur- 
vature effect, as a "straightened out" box, with an orthogonal square cross section of 
size £ (along which the x-y directions lie), and an axial length L (along the z direction) 
embedded in an axial homogeneous uniform magnetic field Bq = Bo e z . This simpli- 
fied geometry allows us to perform simulations with both high numerical resolution and 
long-time duration. 

The dynamics of a plasma embedded in a strong axial magnetic field are well 
described by the equations of reduced MHD (Kadomtsev and Pogutse [13], Strauss 
[14], Montgomery [15]). In this limit the velocity and magnetic fields have only per- 



pendicular components, linked to the velocity and magnetic potentials <p and y by 



u ± = Vx((pe z ), fci = Vx(^ z ). (1) 

Although numerically we advance the equations for the potentials [see 12], in order to 
analyze the linear and nonlinear properties of the system it is convenient to write the 
equivalent equations using the Elsasser variables z ± = u_i±b_i. The more symmetric 
equations, which explicit the underlying physical processes at work, are given in dimen- 
sionless form by: 

£ = _(^.V ± )^ + V,P (2) 

dt ' u ph dz R„ 

dZ ~ (,+ V U- VA dZ ~ JU (- 1 )" +1 y 7 2»„- V p ^ 

-^T — -[Z -V ± Jz - 1 V ± z -V L F (3) 

dt u ph dz R n 

V ± -z ± = (4) 

where P = p + b\/2 is the total pressure, and is linked to the nonlinear terms by 
incompressibility (4): 

2 



P = -L K T )(^)- (5) 



The gradient operator has only components in the x-y plane perpendicular to the axial 
direction z, and the dynamics in the orthogonal planes is coupled to the axial direction 
through the linear terms °< d z . 

We use a computational box with an aspect ratio of 10, which then spans 

0<x,y<l, 0<z<10. (6) 

The linear terms d z are multiplied by the dimensionless parameter va/k p /„ the ratio 
between the Alfven velocity associated with the axial magnetic field v& = Bq/ 'y/4npo, 
and the photo speric velocity u p h- 

Boundary conditions for our numerical simulations are specified imposing the veloc- 
ity potential <p(x,y) in the bottom (z = 0) and top (z = L) planes: 

<p(x,y) = —=^= Y— -^==sm[2n{kx+ly) + 2n^ kl }. (7) 
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These result from the linear combination of large-scale eddies with random amplitudes 
ay and phases %ki (whose values are included between and 1). We excite all the twelve 

independent modes whose wave-numbers are included in the range 3 < [k 2 + 1 2 ) ^ 2 < 4, 
and then normalize the result so that the velocity rms is ~ Umr 1 . 

In terms of the Elsasser variables z ± , to impose a velocity pattern (u*^) at the boundary 
surfaces means to impose the constraint z + +z~ = 2u* ± , and as in terms of characteristics 
(which in this case are simply z ± themselves) we can specify only the incoming wave 



(while the outgoing wave is determined by the dynamics inside the computational box), 
at the top (z = L) and bottom (z = 0) planes the following "reflection" takes place: 



z =—z+2u 
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■z +2u L \ at z = L 
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where u° ± and uh^ are the forcing functions in the respective boundary surfaces. 

At time t = no perturbation is imposed inside the computational box, i.e. b^ = u^ = 
0, and only the axial magnetic field Bo is present: the subsequent dynamics are then the 
effect of the photospheric forcing on the system. 

The linear terms (oc d z ) in equations (2)-(3) give rise to two distinct wave equations 
for the z ± fields, which describe Alfven waves propagating along the axial direction z. 
This wave propagation, which is present during both the linear and nonlinear stages, is 
responsible for the transport of energy at the large perpendicular scales from the bound- 
aries (photosphere) into the loop. The nonlinear terms (z T • V^)z ± are then responsible 
for the transport of this energy from the large scales toward the small scales, where 
energy is finally dissipated, i.e. converted to heat and/or particle acceleration. 

An important feature of the nonlinear terms in equations (2)-(4) is the absence of self- 
coupling, i.e. they only couple counterpropagating waves, and if one of the two fields 
z ± were zero, there would be no nonlinear dynamics at all. This is at the basis of the 
so-called Alfven effect (Iroshnikov [16], Kraichnan [17]), that ultimately renders the 
nonlinear timescales longer and slows down the dynamics. 

From this analysis it is clear that three different timescales are present: Xa, x p h and x n \. 
Xa = L/va is the crossing time of the Alfven waves along the axial direction z, i.e. the 
time it takes for an Alfven wave to cover the loop length L. x p ) x ~ 5 m is the characteristic 
time associated with photospheric motions, while x n \ is the nonlinear timescale. 

For a typical coronal loop Xa <C x p h, and for this reason we consider a forcing which 
is constant in time, i.e. for which formally x p ) % = °°. 

In the RMHD ordering the nonlinear timescale x n \ is bigger than the Alfven crossing 
time Xa - This ordering is confirmed and maintained troughout our numerical simulations. 

The length of a coronal section is taken as the unitary length, but as we excite all the 
wavenumbers between 3 and 4, and the typical convection cell scale is ~ 1000 km, this 
implies that each side of our section is roughly 4000 km long. Our grid for the cross- 
sections has 512x512 grid points, corresponding to ~ 128 2 points per convective cell, 
and hence a linear resolution of ~ 8 km. 

Between the top and bottom plate a uniform magnetic field B = B§e z is present. The 
subsequent evolution is due to the shuffling of the footpoints of the magnetic field lines 
by the photospheric forcing. 

In this section we present the results of a simulation performed with a numerical 
grid with 512x512x200 points, hyper-diffusion (n = 4) with R4 = 10 19 , and the Alfven 
velocity va = 200 km s^ 1 corresponding to a ratio v^/w^ = 200. The total duration is 
roughly 500 axial Alfven crossing times (Xa = L/va). 

Plots of the total magnetic and kinetic energies 
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FIGURE 1. High-resolution simulation with v A /uph = 200, 512x512x200 grid points and R4 = 10 19 . 
Left: Magnetic {Em) and kinetic (E K ) energies as a function of time (t a = L/v A is the axial Alfvenic 
crossing time). Right: The integrated Poynting flux S dynamically balances the total dissipation D. Inset 
shows a magnification of total dissipation and S for 180 < t/x A < 280. 



and of the total magnetic and kinetic dissipation rates 

D M = -^-JdVb ± -V s b ± D K = -^-JdVu ± -V & u ± (11) 

along with the incoming energy rate (Poynting flux) S, are shown in Figure 1. At the 
beginning the system has a linear behavior [see 12], characterized by a time linear 
growth rate for the magnetic energy, the Poynting flux and the electric current, until time 
t ~ 6 Xa, when nonlinearity sets in. The magnetic energy is bigger than the kinetic energy, 
this is the natural result of the field line bending due to the photospheric motions both 
in the linear and nonlinear stages. More formally this is a consequence of the fact that, 
while on the perpendicular magnetic field no boundary condition is imposed, the velocity 
field must approach the imposed boundary values at the photosphere both during the 
linear and nonlinear stages. 

After this time, in the fully nonlinear stage, a statistically steady state is reached, in 
which the Poynting flux, i.e. the energy that is entering the system for unitary time, 
balances on time average the total dissipation rate {Dm + Dk). As a result there is no 
average accumulation of energy in the box, beyond what has been accumulated during 
the linear stage, and a detailed examination of the dissipation time series (see inset in 
Figure 1) shows that the Poynting flux and total dissipations are decorrelated around 
dissipation peaks. 

Figure 2 shows the energy spectra. The spectral index for total energy fits well the 
—2 value. We have shown [see 12] that this spectral index strongly depends on the 
ratio VA/uph, i.e. on the relative strength of the axial magnetic field. At lower values 
correspond flatter spectra, with an index close to —5/3, while to higher values of the 
magnetic field the spectra steepens up to ~ —5/2 for va/w p /, ~ 1000. 
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FIGURE 2. Left: Magnetic, kinetic and total energy spectra averaged in time over <~ 500 T^. The total 
energy spectrum fits a kj 2 power law. Right: Ratio between cross helicity H c and total energy E as a 
function of time. 



CONCLUSION AND DISCUSSION 

The fact that at the large orthogonal scales the Alfven crossing time Ta is the fastest 
timescale, and in particular it is smaller than the nonlinear timescale T M / (which can be 
identified with the energy transfer time at the driving scale), implies that the Alfven 
waves that continuously propagate and reflect from the boundaries toward the interior 
are basically equivalent to an anisotropic magnetic forcing function that stirs the fluid, 
whose orthogonal length is that of the convective cells (~ 1000£:m) and whose axial 
length is given by the loop length L. 

Recently a lot of progress has been made in the understanding of turbulence for 
an MHD system embedded in a strong magnetic field (Ng and Bhattacharjee [18], 
Goldreich and Sridhar [19], Sridhar and Goldreich [20]). As a coronal loop is threaded 
by a strong magnetic field, it is no surprise that the nonlinear dynamics is described by 
weak MHD turbulence. 

The spectra that we have found can be easily derived by order of magnitude consid- 
erations. A characteristic of anisotropic MHD turbulence is that the cascade takes place 
mainly in the plane orthogonal to the DC magnetic guide field (Shebalin et al. [21]). Di- 
mensionally, and integrating over the whole box, the energy cascade rate may be written 
as 

e^Lp^s (12) 
n 

where Szx is the rms value of the Elsasser fields z ± = u\_ ± b\_ at the perpendicular 
scale A. Given the simmetry of the system it is expected and confirmed numerically (see 
Figure 2) that cross helicity is zero, hence bz\ ~ 8z^ ~ 8zx . p is the average density 
and Tx is the energy transfer time at the scale A, which is greater than the eddy turnover 
time %x ~ A/ 8zx because of the Alfven effect [16, 17]. 
In the classical IK case, 

7X~T;A (13) 



More generally, however, as the Alfven speed is increased nonlinear interactions become 
weaker. Simply from dimensional considerations as the ratio %x j %a is dimensionless and 
smaller than 1, we can suppose that the energy transfer time scales as 



rA ~ TA 0l) ' with (i4) 

where a is the scaling index (note that a = 1 corresponds to standard hydrodynamic 
turbulence). 

The energy transfer rate (12) is then given by 

Considering the injection scale A ~ £±, eq. (15) becomes 

On the other hand the energy injection rate is given by the Poynting flux integrated across 
the photospheric boundaries: £,-„ = pVAfdau p h -b±. Considering that this integral is 
dominated by energy at the large scales, due to the characteristics of the forcing function, 
we can approximate it with 

£in ~ P £ 2 1 v A u ph 8ze ± , (17) 

where the large scale component of the magnetic field can be replaced with dzi L because 
the system is magnetically dominated. 

The last two equations show that the system is self-organized because both £ and £ ; „ 
depend on 8ze ± , the rms values of the fields z^ 1 at the scale t\_; the internal dynamics 
depends on the injection of energy and the injection of energy itself depends on the 
internal dynamics via the boundary forcing. 

In a stationary cascade the injection rate (17) is equal to the transport rate (16). 
Equating the two yields for the amplitude at the scale 



8z 
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Uph \ Lu ph J 

Substituting this value in (16) or (17) we obtain for the energy flux 



(18) 



£ *~^ B H^J • (19) 

where va = Bq/ ^/Anp. This is also the dissipation rate, and hence the coronal heating 
scaling. A dimensional analysis of eqs. (2)-(4) reveals [see 12] that the only free pa- 
rameter is / = £^VA/Lu p h, so that the scaling index a (14), upon which the strength 
of the stationary turbulent regime depends, must be a function of / itself, and we have 
determined its value computationally [12]. 



Identifying, as usual, the eddy energy with the band-integrated Fourier spectrum 
8z^ ~ k^Ek ± , where k\_ ~ from eq. (15) we obtain 

3 a +2 

E k± ock ± ^, (20) 

where for a = 1 the —5/3 slope for the "anisotropic Kolmogorov" spectrum is recov- 
ered, and for a = 2 the —2 slope. At higher values of a correspond steeper spectral 
slopes up to the asymptotic value of —3. 

It has been shown computationally and analytically (Ng and Bhattacharjee [18], 
Galtier et al. [22]) that the scattering of Alfven waves with random amplitudes in the 
weak MHD turbulence regime gives rise to the k^ 2 spectrum. Our MHD simulations 
differ from [18] as we integrate forward in time the reduced MHD equations, and the 
system is not three-periodic. While we confirm the presence of the kj 2 spectrum, steeper 
spectra are found up to ~ k^ 3 , and they are clearly linked to the strength of the axial 
field Bq. Future analytical and computational investigation might clarify the physical 
origin for the steeper spectra, their relation to boundary conditions and lying-tying, and 
possibly give an explicit analytical formula for a(f). 
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